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Abstract — The range of possibilities to investigate the 
biomechanical behavior of human biological structures is 
manifold. The classical investigation methods such as e.g. 
experimental studies and image processing are complemented 
by the biomechanical computer modeling. These different 
investigation methods represent their own fields of research, 
but a comparison of the results can be used to validate the 
applied method. Also explored material properties of different 
structures can be used, for example in the computer modeling, 
as input parameters. In order to ensure a realistic modeling, the 
knowledge of the influence of the structure-specific parameters 
on the biomechanical behavior of the entire model system is 
required. Especially the biomechanics of the intervertebral 
discs is seen in the literature as a central component of a spine 
model. Therefore, a detailed analysis of the impact of an input 
parameter variation is necessary. By means of a multibody 
simulation (MBS) model of the human lumbar spine, the impact 
of modified disc stiffness onto the spinal structures has been 
investigated. The used stiffness values are obtained from 
published literature. The model takes the biomechanical 
properties of the spinal structures such as the intervertebral 
discs, the fact joint and the ligamental structures into account 
and has been validated by comparing the results with results 
from appropriate literature. To analyze the effects of different 
input parameters variations on the biomechanical behavior of 
spinal structures, the upright standing is simulated. 

This research project shows that the implementation of 
different input parameters don’t necessarily lead to massive 
changes of the biomechanical behavior of the structures in 
which the input parameter has been varied, but may have a 
greater influence on other modeled structures. 

Index Terms — sensitivity analysis, intervertebral disc, MBS 
model, lumbar spine. 


I. INTRODUCTION 

There are different approaches to identify spinal disorders: 
For example, medical imaging, experimental investigations 
or computer modeling. The most common medical imaging 
procedures are MRI and CT- imaging [1], [2], [3], [4]. With 
the help of these imaging techniques damages of soft tissue 
can be diagnosed and the bone constitution can be assessed. 
Another approach to identify spinal disorders is the 
experimental investigation. In in vivo studies a pressure 
sensor is inserted in a selected disc of a living person, so that 
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the intradiscal pressure can be measured [5], [6], [7], [8]. In in 
vitro studies the kinematic quantities influenced by external 
forces or torques of autopsy can be determined [9], [10], [11]. 
A further method to identify spinal disorders is computer 
modelling. It is distinguished between the multibody 
simulation (MBS) modeling and finite element (FE) 
modeling. Each body of a FE model is divided into smaller 
sub-units. For each element the node displacement and the 
related change in tension is calculated by taking the specific 
material laws into account [12]. However, one disadvantage 
of FE modeling is the required computational time of the 
system, which is relatively long. According to Berkley [13] 
the models' accuracy increases with the number of finite 
elements that are used to describe the geometry. But each 
additional element also means an additional computational 
time. Chomphan [14] and Zhang [15] confirm that solving 
large numbers of FE equations leads to an enormously 
time-consuming calculation. The accuracy of the system and 
the expected computing time must therefore be carefully 
matched. A much faster method is the MBS. In MBS 
modeling the bony structures of the vertebral bodies are 
assumed to be rigid and thus not deformable. The individual 
bodies are linked through massless joints or force elements 
[16]. The acting external force activates the kinematics of the 
model, which is defined as a system of coupled differential 
equations. By numerical integration the kinematic variables, 
the transmitted forces and the torques are calculated. 

One difficulty in modeling is that the input parameters, 
which are partly based on data from already published 
literature, may differ from one to the next publication. For 
example, the values for the stiffness of the intervertebral 
discs are dissimilar in the publication of Lavaste [17]. 

Due to the short computation time of MBS modeling, it is 
possible to analyze the impact of different input parameters 
on a broad spectrum of possible variations. This study 
investigates in which way the results depend on the input data 
that are used to define the material properties of the model. 
The effects of variations in the intervertebral disc material 
properties and the consequences on the spinal structures are 
analyzed. 

II. Multibody Simulation Model ol the lumbar spine 

A. Surface generation and alignment of the lumbar spine 
model 

The vertebral surfaces of this model are based on 
CT-images of artificial vertebrae, whose size correspond to 
the average size of the vertebrae of Europeans. Plugins are 
developed to segment and to visualize the data sets (Fig. 1) 
and make them available for simulation [18], [19]. 
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Fig. 1 Example for a segmented and visualized data set 


using self-developed plugins 

The vertebral bodies are arranged so that the spinal 
alignment fulfills the characteristics of a well-balanced spine 
according to Roussouly [20]. All the conditions of a 
well-balanced spine are fulfilled in this model: the sacral 
slope is between 35° und 45°, the apex of the lumbar lordosis 
is located at the center of the L4 vertebral body, the lower 
lordosis is more prominent, the inflexion point is located at 
the thoracolumbar junction, an average of four vertebral 
bodies constitutes the arc of lordosis, and the average global 
lordosis angle is 61°. 

B. Spinal structures: biomechanical properties 

The lumbar MBS model (Fig. 2) consists of vertebrae L1-L5, 
the os sacrum, and the os ilium. The rigid bodies are 
connected by joints that are located in the middle between 
two vertebral bodies. At this point the forces and torques can 
be transmitted. The disc force is calculated by an equation, 
which is composed of a geometry-based stiffness and a 
damping term (1). The geometry-based stiffness term is 
composed of the stiffness c, the cross-sectional area CSA as a 
unitless factor and the deformation of the intervertebral discs 
Ar. The unitless factor CSA is included in the stiffness term to 
take the effects of the different disc sizes onto the disc 
properties into account. The damping term depends on the 
damping d and the velocity Ar’. 


F = c • CSA • Ar + d • Ar’ (I) 


The transmission of torques is based on experimentally 
determined curves for all three motion axes [21]. 



Fig. 2 MBS-model of the lumbar spine 


Furthermore, the ligaments lig. longitudinale posterius 
(PLL), lig. longitudinale anterius (ALL), lig. flavum (LF) and 
lig. interspinale (ISL) as well as the lig. supraspinale (SSL) 
and the lig. intertransversarium (ITL) and ligg. iliolumbale 
are implemented in the model. A ligament is spanned 
between to marker points. As a ligament can only be defined 
between two points in the simulation, broad ligament 
structures are realized by a bundle of several fibres. The 
mechanical behaviour of the ligaments is also based on 
characteristic curves, which describe the 
force-deformation-relation of the individual ligaments [22]. 

In addition ten facet joints are included with a contact 
modeling. If the two corresponding facet joint surfaces are in 
contact, a force is developed in opposite direction of the 
movement (2). 


'Fy' 


c y • Ar + d v • Ar‘: c y <0; Ar<0; Ar‘< 0 

c 


c y -Ar : c v <0; Ar<0; Ar‘> 0 

Fy 

= < 

0 : c v < 0; Ar > 0; 

F, 


0 

UJ 


0 

^ > 


A more detailed description of the biomechanical 
properties of the structure can be taken from [26]. 

C. Validation of the model 

The model validation was performed by comparing the 
simulation results with FE results and in vivo data from the 
literature [23], [24], [25], [26]. As an example, the pressure of 
the intervertebral discs of the each functional spine unit 
(FSU) is shown in Fig. 3. 
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Fig. 3 Intervertebral disc pressure in comparison 
Comparing the intervertebral pressure of the FSUs, it can 
be seen that the values are not exactly in the same order 
concerning magnitude. Particularly obvious are the 
differences in the results for pressure in the FSUs L5-Sac via 
MBS modeling (Bauer) and the finite element modeling 
(Rohlmann). 

The rotational behavior of the discs and the loads of the 
facet joints are also validated with corresponding results from 
literature. A detailed description of the validation process and 
relevant conclusions about possible causes of discrepancies 
in results are shown in [22]. 


III. Realization of the disc parameter variation 

As already mentioned, different biomechanical parameters 
for the intervertebral disc stiffness can be found in literature. 
After the implementation of the different values, its effects on 
the spinal structures are examined. Lavaste [17] specifies his 

o 

experimental obtained stiffness value c with c=8*10 N/m, 

o 

his modeled stiffness value with c=9*10 N/m and indicates 
further stiffness values from experimental studies of other 
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scientists with c=10*10 8 N/m (stiff_Lavaste_Markolf), 
c=13*10 8 N/m (stiff_Lavaste_Panjabi), c=5.2*10 8 N/m 
(stiff_Lavaste_T encer). 

To analyze the changes in the biomechanical behavior of 
the spinal structures, two series of studies are performed. In 
the first for all functional spinal units the same stiffness 
parameters are used. In the second study all possible 
combinations of the different stiffness parameters and all 
functional spinal units constitute the basic parameter 
configuration. With a total number of five different stiffness 
input parameters and five functional spinal units a possibility 
of 3125 combinations arises. For this number of 
combinations the loads of the intervertebral discs are 
calculated. 

The most natural case, the upright position, was simulated 
as a load case. This means that an external force of 500N, 
which corresponds to the weight of the upper body, is applied 
in vertical direction on the top of the surface of vertebra LI. 
By this external force, the spinal structures are brought out of 
their equilibrium state before being transferred to a new one. 
This new equilibrium state is considered in the following 
results. 


IV. Results 

A. 1 st study: same stiffness parameters for all FSUs 

Considering the loads of the intervertebral discs, it can be 
seen that the forces in different FSUs are in a range between 
400N and 500N (Fig. 4). In comparison the intervertebral 
discs of FSUs L5-Sac are loaded the lowest. The disc forces 
increase with in the subsequent FSUs. The largest force is 
developed in the intervertebral disc L2-L1. A possible reason 
could be the alignment of the vertebral bodies L2 and LI and 
thus, the direction orientation of the force vector of the 
intervertebral disc L2-L1. In general, the smaller the 
inclination value of the intervertebral disc the higher the 
vertical component and the lower the horizontal component 
force of the intervertebral disc. In the presented MBS-model 
particularly the intervertebral disc of the FSU L2-L1 is 
slightly inclined so that the horizontal component of the 
intervertebral disc force is relatively small and the vertical 
force component correspondingly higher. A detailed 
explanation concerning the relationship between the 
alignment of the vertebrae and the force components can be 
seen in [22]. Within a FSU the reaction forces of the 
intervertebral disc are almost identical. Only small variation 
in disc forces can be reported in the FSUs L4-L5-L2-L3. The 
deviations are for this FSUs under 3%. 


□ stiffLavasteTencer ■ stifFLavasteexp 
8 stiff Lavaste model CD stiff Lavaste Markolf 



Fig. 4 Intervertebral disc force in comparison 


The influence of different stiffness parameters on the 
extension and flexion movement of the intervertebral discs is 
nearly similar. Figure 5 shows on the x-axis the extension and 
flexion movement in degrees and on the y-axis the FSUs. The 
positive x-values describe the ventral directed flexion and 
negative x-values the dorsal directed extension movement. 
Within the individual FSUs the rotational directions are 
identical. But it should be noted that the simulated loading 
case, the upright standing position, generally causes only 
very small rotations. This is due to the fact that this modelled 
lumbar spine is "well-balanced” and thus, the alignment of 
the vertebrae is physiologically optimally. The comparison of 
the amounts of rotation of the individual FSUs with different 
stiffness parameters shows that the smaller stiffness value 
stiff_Lavaste_Tencer (c=5.2*10 8 N/m) causes smaller 
rotation in the FSUs L5-Sac, L3-L4 and L2-L3 than using the 

o 

higher stiffness value stiff_Lavaste_Panjabi (c=13*10 N/m). 
The situation of the FSUs L4-L5 is exactly reversed. The 
higher stiffness parameters of Panjabi evoke larger rotations 
in these FSUs, than the stiffness parameters 
stiff_Lavaste_T encer, stiff_Lavaste_exp, 

stiff_Lavaste_model and stiff_Lavaste_Markolf. The largest 
deviations, 26%, can be found for the FSU L2-L3. 
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Fig. 5 Comparison of extension and flexion 
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Fig. 6 shows the loads of the facets joints of the FSU using 
different interdiscal stiffness parameters. In this load case, 
the corresponding facet surfaces touch and according to the 
equation (2) build an appropriate contact force. The reason, 
why only a very small force is built in the FSU L1-L2, is the 
alignment of the facet surfaces with respect to the acting 
external force. In particular, the facet surfaces of the FSU 
L1-L2 are aligned in parallel to the line of action of the 
external force so that the two surfaces slide past each other 
and have only little contact. In this case, a very small contact 
force is build up. It is evident that for all FSUs the facet force 
is higher when using small interdiscal stiffness parameters 
and larger when implementing higher stiffness parameters. 

The stiffness has the biggest impact on the FSU L2-L3. In 
this FSU the smaller stiffness value stiff_Lavaste_Tencer 
also causes higher loads of the facet joints than the higher 
stiffness value stiff_Lavaste_Panjabi. This deviation can be 
justified by the above described characteristic of the 
rotations. For this FSU L2-L3 the extension has a direct 
impact on the loads of the facet joints. Due to higher 
backward tilting the posterior facet joints are more heavily 
loaded in this FSU. 
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Fig. 6 Comparison of the facet loads 


B. 2 nd study: all possible combinations 

The large amount of result data, which corresponds to a 
number of 3125 simulations, allows an even more intensive 
insight into the sensitivity of the model with respect to 
parameter variation. In this process the parameter 
configuration is determined in which the intervertebral discs 
of each FSU are loaded to the maximum (Table 1) and 
minimum (Table 2). In addition, the mean load of the 
intervertebral discs of each FSU and the standard deviation 
are calculated (Fig. 7). 


still strongly influence the load situation of the intervertebral 
discs. 

The mean values of the intervertebral disc force, calculated 
from the 3125 parameter combinations, and the 
corresponding standard deviation are shown in Fig. 7. The 
standard deviation of the intervertebral disc forces is 
calculated as follows 


SD 


Cr \2 


Z i=1 (Xi- x) 
(n) 



with Xi the i-th element of the sample, x the mean value and n 
the numbers of values. 


600 



L5-Sac L4-L5 L3-L4 L2-L3 L1-L2 

Fig. 7 Mean values and standard deviation 


Table 1 Maximum load of the intervertebral discs 



max. F disc 
[N] 

stiffness parameter configuration of the different FSUs 

L5-Sac 

[N/m| 

L4-L5 [N/m] 

L3-L4 [N/m] 

L2-L3 [N/m] 

L1-L2 [N/m] 

L5-Sac 

435,6 

10- 10* 

9- I0 8 

8- 10 8 

DC 

o 

o 

9- 10 8 

L4-L5 

462,9 

10- 10 8 

9 10 8 

oc 

O 

OO 

OO 

o 

o 

9 - 10 8 

L3-L4 

495,0 

10- 10 8 

13 • 10 8 

10 - 1 o 8 

10- 10 8 

5,2- 10 8 

L2-L3 

487,6 

10' 10 8 

9- 10 8 

8 • 10 8 

10- 10 8 

9- 10 8 

L1-L2 

499,7 

9- 10 8 

8- 10 8 

10- I0 8 

8 • 10 8 

10- 10 8 


Comparing the maximum loads of the intervertebral discs 
with those who were reached by the parameter configuration 
in Fig. 4, it can be seen that the load of the FSUs are nearly 
equal. 

Table 2 Minimum load of the intervertebral discs 



min. F UiiC [N] 

stiffness parameter configuration of the different FSUs 

— r 

z V 
S' & 

3 CS 

— n 

L4-L5 [N/m] 

L3-L4 [N/m] 

L2-L3 [N/m] 

L1-L2 [N/m] 

L5-Sac 

266,0 

9- 10 8 

13 • 10' 

13 • I0 8 

5,2- 10 8 

10* 10 8 

L4-L5 

299,6 

9- 10 8 

8- 10 8 

13 • 10 8 

LA 

N> 

O 

OO 

5,2- 10 8 

L3-L4 

345,7 

9- 10 8 

8- 10 8 

13 • 10 8 

5,2 • 10 8 

5,2- 10 8 

L2-L3 

335,8 

10- 10 8 

8- 10 8 

13 • I0 8 

5,2* 10 8 

5,2- 10 8 

L1-L2 

455,9 

9- 10 8 

13- 10 8 

13 • 10 8 

5,2 • 10 8 

10- 10 8 


By contrast, highly different is the variance between the 
maximum values and the minimum values of Table 1 and 
Table 2. For example, the difference of the disc loading of the 
FSU L5-Sac lies at 169N, which corresponds to a 39% 
increase of the intervertebral disc loading. Conversely, this 
means if selecting the parameter configuration of the case 
“minimum load situation”, small values are calculated for the 
disc forces, but this does not necessarily have to be correct. 
This example shows that, depending on the choice of 
parameters configuration, the input parameter stiffness can 


V. Conclusion 

This study investigating the influence of the intervertebral 
disc stiffness parameter on the biomechanical behavior of the 
lumbar spine shows that some structures are less sensitive to 
changes of the intervertebral disc stiffness than others, such 
as the facet joints. As shown above, within each FSUs the 
facet joints are more than twice as heavily loaded. These 
significant differences in the facet loads, which have been 
purely caused by the variation of the intervertebral discs 
stiffness, underline that the modeling results largely 
dependent of the input parameters. 

Although, according to Fagan [28], the intervertebral disc 
is the most critical component in most finite element models 
of the spine and its representation in the models therefore of 
great importance, it is also essential to evaluate the parameter 
influence of further spinal structures. Therefore, in a further 
study, we will investigate the influence of different 
characteristic curves of the ligaments, which are found in the 
literature. 

In conclusion, the characteristics of the individual spinal 
structures cannot be considered in isolation, but the 
biomechanical behavior of certain structures can influence 
other structures. Currently the presented biomechanical 
behavior of the intervertebral disc is defined by a force law, 
which can be understood as an initial approach. To precise 
the biomechanical properties of the intervertebral disc a 3D 
hybrid model consisting MBS and FE units will be build. 
Lastly, it should be noted that, after further sensitivity 
analyses, we target patient-specific preoperative simulations 
to predict the effects on spinal fusion in overweighed and 
obese patients and to identify the best possible surgical 
option. 
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